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Abstract 

We consider an infinite time horizon spatially distributed optimal harvesting prob¬ 
lem for a vegetation and soil water reaction diffusion system, with rainfall as the main 
external parameter. By Pontryagin’s maximum principle we derive the associated four 
component canonical system, and numerically analyze this and hence the optimal control 
problem in two steps. First we numerically compute a rather rich bifurcation structure 
of flat (spatially homogeneous) and patterned canonical steady states (FCSS and PCSS, 
respectively), in ID and 2D. Then we compute time dependent solutions of the canonical 
system that connect to some FCSS or PCSS. The method is efficient in dealing with 
non-unique canonical steady states, and thus also with multiple local maxima of the ob¬ 
jective function. It turns out that over wide parameter regimes the FCSS, i.e., spatially 
uniform harvesting, are not optimal. Instead, controlling the system to a PCSS yields a 
higher profit. Moreover, compared to (a simple model of) private optimization, the social 
control gives a higher yield, and vegetation survives for much lower rainfall. In addition, 
the computation of the optimal (social) control gives an optimal tax to incorporate into 
the private optimization. 

Keywords: distributed optimal control, bioeconomics, optimal harvesting 
MSC: 49J20, 49N90, 35B32 


1 Introduction and main results 


Vegetation patterns such as spots and stripes appear in ecosystems all over the world, in 
particular in so called semi arid areas DBC^OS] . Semi arid here means that there is enough 
water to support some vegetation, but not enough water for a dense homogeneous vegetation. 
Such systems are often modeled in the form of two or more component reaction-diffusion sys¬ 
tems for plant and water densities, with rainfall R as the main bifurcation parameter, and the 
patterns are attributed to a positive feedback loop between plant density and water infiltra¬ 
tion. Starting with a homogeneous equilibrium of high plant density for large i?, stationary 
spatial patterns appear as R is lowered, often following a universal sequence |GRS14| . This 
may lead to catastrophic, sometimes irreversible, regime shifts, where a patterned vegetation 
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suddenly dies out completely, leaving a desert behind as R drops below a certain threshold, 
or as some other parameter such as harvesting or grazing by herbivores is varied. There 
is a rather large number of specific models, each agreeing with field observations in vari¬ 
ous parameter regimes, see e.g., |SZvHMdl] iHRvdB^Ol , and |RDdRvdK04| for a review, or 
|SBB^09[ IZKY^IS) and the references therein for further reviews including more recent work 
on early warning signals for desertification and other critical transitions, usually associated 
with subcritical bifurcations. 

A so far much less studied problem is the spatially distributed dynamic optimal control 
(OC) of vegetation systems by choosing harvesting or grazing by herbivores in such a space 
and time dependent way that some economic objective function is maximized. Following 
|BX10| we consider an infinite time horizon OC problem for a reaction-diffusion system for 
vegetation biomass and soil water, which is roughly based on |HRvdB^0l| . Related optimal 
control problems have also been considered in |BX08[ |XeplO[ IBEGX13) , with the focus on 


the so called “Optimal Diffusion Instability” of flat canonical steady states (FCSS, see [1.1 


for OC related definitions), which similar to a Turing bifurcation yields the bifurcation of 
spatially patterned canonical steady states (PCSS) from FCSS. However, in these works, and 
in jBxini, only few PCSS have been actually calculated numerically, and no canonical paths, 
i.e., time dependent solutions of the canonical system. Finite time horizon cases recently 
have been considered in |CPB12[ IACKLTT^ IADSI^ . mostly focusing on theoretical aspects, 
and on problems with control constraints, which altogether gives a rather different setting 


than ours; see also Remark 1.3 


Here we apply the numerical framework from |GU15[ lUeclSj to the so called canonical 
system for the states and co-states, derived via Pontryagin’s Maximum Principle. First we 
calculate (branches) of canonical steady states (CSS), including branches of PCSS, and in a 
second step canonical paths connecting to some CSS. Our main result is that in wide ranges of 
parameters, in particular for low rainfall i?, FCSS and their canonical paths are not optimal, 
and that controlling the system to a PCSS yield a higher profit than uniform harvesting. 
Thus, this seems to be the first example of a bio-economically motivated optimal control 
problem, where the global bifurcation structure of CSS has been computed in some detail, 
showing multiple CSS at fixed parameter values and with dominant non-spatially homoge¬ 
neous steady states, and where moreover canonical paths to such PCSS have been computed, 
with significant gains in welfare. We also compare these results to results for the same system 
with (initially) no external control, i.e., the system with private optimization, and find that 
the optimal control significantly increases the profit, and supports vegetation at significantly 
lower rainfall levels, which again has important welfare implications. Remarkably, in our 
system the co-state of the vegetation can be identified with a tax for the private optimiza¬ 
tion problem, and thus, by solving the canonical system we find an optimal space and time 
dependent tax for the private optimization problem. 

A standard reference on ecological economics or “Bioeconomics” is |Cla90| . including a 
very readable account, and applications, of Pontryagin’s Maximum Principle in the context 
of ODE models; see also |LW07j . A review of management rules for semi arid grazing systems 
including comparison to real data, and making plain the importance of such rules, is given in 
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However, the (family of) models discussed there consist of discrete time evolutions 
without spatial dependence, and with focus on the rainfall i? as a time-dependent stochastic 
parameter; here we consider a deterministic PDE model. Thus comparison between the two 
(classes of) models is difficult, but it should be interesting to include spatial dependence into 
the model in |QB12 . 


Finally, in, e.g., |Neu03j . |DL09| . stationary spatial OC problems for a fishery model are 
considered, including numerical simulations, which correspond to our calculation of canonical 
steady states for our model © below. The results of |Neu03[[DL09| show that for their models 
it is economically optimal to provide “no-take” marine reserves. This is similar to our finding 
of optimal patterned canonical steady states, but here we go beyond the steady case with 
the computation of optimal paths. This in particular tells us how to dynamically control the 
system to an (at least locally) optimal steady state. Moreover, even for the steady states there 
are important technical differences between the OC problems considered in |Neu03[ IDL09| 
with unique positive canonical steady states, and the OC problems considered here, which 
in relevant parameter regimes are distinguished by having many canonical steady states, of 


which several can be locally stable (see Remark 1.2 for definitions, and, again. Remark |l.3[ ) 
Thus, management rules for our system are considerably more complicated than those in, 
e.g., |Neu03[ rDL09| . as they must take the different CSS into account, and given an initial 
state, must decide to the domain of attraction of which CSS it belongs. Partly to illustrate 
this point, we also compute a Skiba (or indifference) point |Ski78j between different CSS in 


^2.3.2 


In the remainder of this Introduction we explain the model and the use of Pontryagin’s 
Maximum Principle to derive the canonical system (0), explain the basic idea of the nu¬ 
merical method (M.2), and summarize the results ( ^1.3[ ). In ^we present the quantitative 
results, and in ^we give a brief further discussion. The method has been implemented in 
Matlab as an “add-on” package p2pOC to the continuation and bifurcation software pde2path 
|UWR14] . and the matlab functions and scripts to run (most of) the simulations, the under¬ 
lying libraries, manuals of the software, and some more demos can be downloaded at |Uec H- 

Acknowledgment. I thank D. Grafi, ORCOS Wien, for valuable comments on the economic 
terminology, and the anonymous referees for valuable questions and comments on the first 
version of the manuscript. 

1.1 Problem setup 

In dimensionless variables the model for vegetation and soil water with social optimization 
from |BXI0j is as follows. Let C be a one-dimensional (ID) or two-dimensional 
(2D) bounded domain, v = v{x^t) the vegetation density at time t > 0 and space x G D, 
w = w{x^t) the soil water saturation, 

E = E{x^t) the harvesting effort (control), H{v^E) = v^E^~^ the harvest, and 


I r 

= Pi 


^-Pt 


/ Jc{v{x^t)^E{x^t)) dxdt 

Jn 


( 1 ) 
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the (spatially averaged) objective function, where p > 0 is the discount rate, and 

Jc{v,E)=pH{v,E)-cE (2) 

is the (local) current value profit. This profit depends on the price p, the costs c for harvesting, 
and E in a, classical Cobb-Douglas form with elasticity parameter 0 < cr < 1. The problem 
reads 


V{vo^wo) = maxJ{vo^wo^E)^ where (3a) 

dtv = diAv + [gwv^ — d{l + 6v)]v — H{v^ E)^ (3b) 

dtw = d2Aw + i?(/3 + ^v) - {tuV + ruj)w, (3c) 

{v,w)\t=o = ('?;o,'^^o)- (3d) 


The parameters and default values are explained in Table together with brief comments 
on the terms appearing in Q. One essential feature of ([^,c) is the positive feedback loop 
between vegetation v and water w. Clearly, the vegetation growth rate gwv^ increases with 
re, but the vegetation v also has a positive effect on water infiltration, for instance due to 
loosened soil, modeled by R^v in 

Remark 1.1. In |BX10j . E is also refered to as the grazing by herbivores, and Q is called a 
semi arid grazing system. This seems somewhat oversimplified, because from a practical point 
of view, controlling in detail the movement and grazing behavior of, e.g., cattle, certainly is 
a more complicated problem than controlling genuine harvesting. Thus, henceforth we stick 
to calling E the harvesting effort. J 


We complement © with homogeneous Neumann boundary conditions (BC) — 0 
and djyW = 0 on where u is the outer normal. The discounted time integral in 0 is 
typical for economic (here bioeconomic) problems, where “profits now” weight more than 
mid or far future profits. More specifically, p corresponds to a long-term investment rate. We 
normalize Jca by |^| easier comparison between different domains and space dimensions. 
In (bio)economics, the control chosen externally by a “social planner” to maximize the 
social value J, is often called social control, as opposed to private optimization, see (19) 
below. Finally, the max in ([^) runs over all admissible controls E] essentially this means that 
E G L^([0, oc) X yi,M), where moreover implicitly we have the control constraint E > 0^ and 
state constraint v^w >0 for the associated solutions of ([^,c). However, in our simulations 


these constraints will always naturally be fulfilled, i.e., inactive^ see also Remark 1.3 


Introducing the costates (A,/i) = (A,/i)(x,t) and the (local current value) Hamiltonian 


T-L{v^ le. A, /i, E) — Jc{v^ E) -\- X \diAv + {gwv^ — d{l + 6v))v — i7] 

+ fj, [d2Aw + R{l3 + ^v) - (ruV + ryj)w], 

by Pontryagin’s Maximum Principle for R = df with the spatial integral 

T-Lit) = / 'H{v{x,t),w{x,t),X{x,t),iJ,{x,t),E(x,t))dx, 

Jn 


( 4 ) 

( 5 ) 
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param. 

meaning 

default values 


coefficient and exponent in plant growth rate gwv^ 

g = 0.001, T/ = 0.5 

d,S 

coefficients in plant death rate d{l + 6v) 


d = 0.03, (5 = 0.005 


coefficients in the infiltration function ^ 

/? = 0.9,,^ = 0.001 

R 

rainfall parameter, used as main bifurcation parameter 

between 4 and 100 

ru,rw 

water uptake and evaporation parameters in the water 

= 0.01, = 0.1 


loss rate TyV + Tyj 




di,2 

diffusion constants for vegetation and water (resp.) 

di = 0.05, d2 = 10 

P 

discount rate 



p = 0.03 

c,p, a 

(economic) param. in the harvesting H{v 

, E) = v^E^ ^ 

c = l,p = 1.1, a = 0.3 


and in the value Jc{v^ E)—pH[v^ E)—cE 




(p and p used as bifurcation param. in ^ 

2.5 

) 



Table 1: Dimensionless parameters and default values in Q; see |BX10j for further comments 
on the modeling. In particular, following |BX10| §4.2] we have a rather larger d 2 . 


an optimal solution A,/i) has to solve the canonical system (CS) 

dtv — d\T-L — di/S.v + [gwv^ — d{l + 6v)]v — i7, (6a) 

dfW = d^T-L = d2Aw + i?(/3 + - {tuV + (6b) 

dfX — p\ — dyT-L — p\ — pav^~^ — A \g{ri + l)wv'^ — 2d6v — d — E^~^] (6c) 

- - Vuw) - diAA, 

dtp = pp- dyjT-L = pp- + p{ruV + r^j) - d2Ap, (6d) 


where E = argmax^ A, /i, E)^ which is obtained from solving OeR = 0 for E^ giving 

«>a) 

The costates (A, p) also fulfill zero flux BC, and derivatives like dyT-L etc are taken variationally, 
i.e., for T-L. For instance, for $('?;, A) = XAv we have $('?;, A) = f^XAvdx = f^(AX)vdx by 
GauB’ theorem, hence 6y^{v^ X)[h] = f (AX)hdx, and by the Riesz representation theorem 
we identify Sy^(v, A) and hence dy^(v, A) with the multiplier AA. Moreover, we used the so 
called intertemporal transversality conditions 

lim e“^^ / v(x,t)X(x,t)dx = Oand lim e“^^ / w(x,t)p(x,t)dx = 0, (6f) 

t^OO Jq t^OO Jq 

which is justified since we are only interested in bounded solutions. 

For convenience setting 

u := ('u, w^X^ p) : Q X [0, oo) ^ M^, (7) 
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we collect & e) and the boundary conditions into 




^ di 

0 

0 

0 ^ 


dtu = —G{u) \— 

VAu + f{u), V = 

0 

d2 

0 

0 

, (8a) 

0 

0 

—di 

0 



(o 

0 

0 

-^ 2 ) 


djyU = 0 on 

{v,w)\t=Q = 





(8b) 


In (8b) we only have initial conditions for the states, i.e., half the variables, and (8a) is 
ill-posed as an initial value problem due to the backward diffusion in (A, /i). Thus, below we 
shall further restrict the transversality condition ([^ to requiring that u{t) converges to a 
steady state, i.e. a solution of 


G{u) = 0, diyU = 0 on dQ. (9) 

Remark 1.2. (Definitions and notations) A solution u of the canonical system ^ is 
called a canonical path^ and a solution u of (§ is called a canonical steady state (CSS). 
With a slight abuse of notation we also call E) with E given by #) a canonical path, 
suppressing the associated co-states A, /i. In particular, if is a CSS, so is (i), w^E). A CSS u 
is called flat if it is spatially homogeneous, and patterned otherwise, and we use the acronyms 
FCSS and PCSS, respectively. For convenience, given t u{t) we also write 


J(7i) := J{vq^wq,E), 


( 10 ) 


with ('L’ 0 ,'^ 0 ) and E (via ^)) taken from u. A canonical path u (or (v^w^E)) is called 
optimal if there is no canonical path starting at the state values and yielding a 

higher J than J{u). As a special case, a CSS u = is called optimal if there is 

no canonical path starting at (v^w) and yielding a higher J than J{v^w^E). We use the 
acronyms OSS for any optimal CSS, and FOSS and POSS for optimal ffat or patterned CSS 
respectively. An OSS u is called locally stable if for all admissible ('L’o^'^o) close to {v^w) 
there is an optimal path u with and lim^^ooSimilarly, u is called globally stable 
if for all admissible ('L’o^'^o) fhe associated optimal path has lim^^ooFinally, the 
domain of attraction of a locally (or globally) stable OSS u is defined as the set of all (vq^ wq) 
such that the associated optimal path u fulfills lim^^oo'^(^)='^- See also |GU15| for more 
formal definitions, and further comments on the notion of optimal system, and, e.g., the 
transversality condition (6f). J 


Remark 1.3. For background on OC in a PDF setting see for instance |Trol0j and the refer¬ 
ences therein, or specifically |RZ99a[ IRZ99b[ ILM01| and [AAClli Chapter5] for Pontryagin’s 
Maximum Principle for OC problems for semilinear parabolic state evolutions. However, 
these works are in a finite time horizon setting, often the objective function is linear in the 

^For “close to” and lim we may use, e.g., the H^{Q) norm, but since all solutions will be smooth, for instance 
as solutions of semilinear elliptic systems with smooth coefficients, we decided to omit the introduction of 
function spaces here. Similarly, (vq^wq) “admissible” should be read as vo(x),wo(x) > 0 for x G P. 
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control, and there are state or control constraints. For instance, denoting the control by 
k = k{x^t)^ often k is chosen from some bounded interval iF, and therefore is not obtained 
from the analogue of ([^), but rather takes values from dK^ which is usually called bang-bang 
control. In, e.g., |CPB12[ IACKLTI3] . some specific models of this type have been studied in 
a rather theoretical way, i.e., the focus is on deriving the canonical system and showing well- 
posedness and the existence of an optimal control. |ADS14| additionally contains numerical 
simulations for a finite time horizon control-constrained OC problem for a three species spa¬ 
tial predator-prey system, again leading to bang-bang type controls. See also |NPS11| and 
the references therein for numerical methods for (finite time horizon) constrained parabolic 
optimal control problems. 

Similarly, the (stationary) fishery problems in |Nen03[[DL09| come with harvesting (i.e. con¬ 
trol) constraints. Moreover, in contrast to our zero-flux boundary conditions ([^) Dirichlet 
boundary conditions are imposed. In |KXL15| the models from |Nen031 IDL09| are extended 
to Robin boundary conditions, and a finite time horizon, with a discounted profit of the form 
J = j^pe~P^h{x^t)u{x^t)dxdt^ where p, p > 0 denote the price and discount rate, h is 
the harvest, and u the (fish) population density, which fulfills a rather general semilinear 
parabolic equation including advection. The first focus is again on well-posedness and the 
first order optimality conditions, and numerical simulations are presented for some specific 
model choices, illustrating the dynamic formation and evolution of marine reserves. However, 
the setting again is quite different from ours, due to the finite time horizon, and since J is 
linear in h and ?i, and since consequently there are constraints on /i, leading to (unique) 
bang-bang controls. 

Here we do not consider (active) control or state constraints, and no terminal time, but 
the infinite time horizon. Our models and method are motivated by |BX08[ IBX10| , which 
also discuss Pontryagin’s Maximum Principle in this setting. We do not aim at theoretical 
results, but rather consider Q after a spatial discretization as a (large) ODE problem, and 
essentially treat this using the notations and ideas from |BPS01| and |GCF^08[ Chapter 7], 
and the aDorithms from |GU15[ lUec m, to numerically compute optimal paths. J 


1.2 Solution steps 

Using the canonical system Q we proceed in two steps: first we compute (branches of) CSS, 
and second we solve the “connecting orbit problem” of computing canonical paths connecting 
to some CSS. Thus we take a broader perspective than aiming at computing just one optimal 
control, given an initial condition (vq^wq)^ which without further information is an ill-posed 
problem anyway. Instead, our method aims to give a somewhat global picture by identifying 
the (in general multiple) optimal CSS and their respective domains of attraction, as follows: 
(i) To calculate CSS (which automatically fulfill ([^) we set up Q as a bifurcation problem 
and use the package pde2path |UWR14( IDRUW14) to find a branch of FCSS, from 
which various branches of PCSS bifurcate. We focus on the rainfall parameter R as 
the main bifurcation parameter, but in ^2^ also briefly discuss the dependence on the 
discount rate p and the price p as bifurcation parameters. By calculating in parallel 
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the (spatially averaged) current value profits 

Jc,a{v, £') := ^ Jc{v{x),E{x)) dx ( 11 ) 

of the CSS we can moreover immediately find which of the CSS maximize Jc^a cind 
hence J — Jc,a/P amongst the CSS. 

(ii) In a second step we calculate canonical paths ending at a CSS (and often starting at 
the state values of a different CSS), and the objective values of the canonical paths. 
Using a Finite Element Method (FEM) discretization, Q is converted into the ODE system 
(with a slight abuse of notation) 




( 12 ) 


where u G is a large vector containing the nodal values oiu — ('u, u;. A, /i) at the n spatial 
discretization points (typical values are n = 30 to n = 100 in ID, and n — 1000 and larger 
in 2D), and M G is the so called mass matrix, which is large but sparse. In ( [l^ 

and the following we mostly suppress the dependence of G on the rainfall R (or the other 
parameters). For (i) we thus need to solve the problem G{u) — 0, which can be considered 
as an elliptic problem after changing the signs in the equations (^ 3 , 4 ( 1 /) = 0 for the costates. 

For (ii) we choose a suitable truncation time T > 0 and replace the transversality condition 
( 6 f) by the condition 


u{T) G Ws{u) and \\u{T) — u\\ small. 


(13) 


where Ws{u) is the stable manifold of the CSS u for the finite dimensional approximation 
(12) of In practice we use 


u(T) G Es{u) and \\u{T) — u\\ small, (14) 

where Es{u) is the stable eigenspace of i.e., the linear approximation of Ws{u) at u. At 
t = 0 we already have the boundary conditions 


(v,w)|t=o = (^’o,^yo) e 


(15) 


for the states. To have a well-defined two point boundary value problem in time we thus need 

dim Es{u) = 2n. (16) 

Since the (generalized, in the sense of M on the left hand side of ([T^) eigenvalues of the 
linearization —duG{u) of ^ around u are always symmetric around p/ 2 , see |GU15[ Appendix 
A], we always have dimEs{u) < 2n. The number 

d{u) — 2n — dimE 5 (?i) (17) 

is called the defect of fi, a CSS u with d{u) > 0 is called defective^ and if d{u) = 0 , then u 
has the so called saddle point property (SPP). Clearly, these are the only CSS such that for 
all ('L’ 05 '^ 0 ) close to (f),u)) we may expect a solution for the connecting orbits problem 

= —G{u)^ {v^w)\t=o = ('L’o^'^o)? G Es{u)^ and \\u{T) — u\\ small. (18) 

See |GU15j for further comments on the significance of the SPP ( [T^ on the discrete level, 
and its (mesh-independent) meaning for the canonical system as a PDE. 









1.3 Results 


The bifurcation diagram (i) for Q turns out to be quite rich, already over small ID domains. 
Thus we mostly focus on ID, but we include a short 2D discussion in §2.4[ Details will be 
given in ^ but in a nutshell we have: In pertinent R regimes there are many CSS, but most 
of them do not fulfill the SPP, and most of these that do fulfill the SPP are not optimal. 
On the other hand, in particular at low R there are locally stable POSS (patterned optimal 
steady states). Before further commenting on this, we briefly review results for the so called 
uncontrolled case. 

In |BX10[ §2] it is shown that the case of private objectives, where economic agents 
(ranchers) are located at each site x, and each one maximizes his or her private profits, leads 
to the system 


dtv = di/S.v + {gwv'^ — d{l + Sv) — A)v, 
dtw = d2Aw + R{f3 + - {tuV + r^)w, 


(19a) 

(19b) 


i.e., the harvest is H = Av^ A > 0 a, fixed parameter. In detail, ranchers with certain property 
rights individually maximizing 7v{v^ E) = — cE leads to 


E — with 7 = ( ^and hence A — ^^ ^. 


( 20 ) 


The same happens in the so called open access case, where agents may harvest freely, giving 

E — with 7 = ( - 1 and hence A = 7 ^“^. On the other hand, (19) can also be seen as 

a “undisturbed Nature” case with modified vegetation death rates d = d + A and S = 

For the economic parameters (c,p, a) = (1,1.1, 0.3) from Table[^we obtain A = 0.543, which 
is rather large compared to the original d = 0.03 from Table 

The bifurcation picture for steady states of (19) is roughly similar to that of ([^, but there 
are also significant differences, and we altogether summarize our results as follows: 

(a) For large R there is a “high vegetation” FCSS, which is a globally stable FOSS (flat 
optimal steady state). 

(b) For smaller R the FCSS from (a) loses optimality, and there bifurcate branches of (lo¬ 
cally stable) POSS (Patterned Optimal Steady States). In particular, we can calculate 
canonical paths from the state values of the FCSS to some POSS which increase the 


profit (up to 40%, see 32). 


(c) The uncontrolled flat steady states (FSS) of (19) only exist for much larger R values 
than the FCSS. At equal R, the profit J (or equivalently the discounted value J/p) of 
the FSS is much lower (e.g., one tenth, see Table bottom center) than the value of 
the FCSS of ([§. 


(d) For the initial value problem (19) we may consider the stability of steady states, while 
CSS at best have the SPP. It turns out that the FSS branch loses stability at a much 


larger value of i? to a patterned steady state of (19), than the FOSS loses optimality in 


the optimal control problem (©• Thus, in the uncontrolled problem pattern formation 
sets in at larger R. 
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Roughly speaking, (b) means that at low R it is advantageous to restrict harvesting to cer¬ 
tain areas. This is similar to the marine reserves in the fishery models in |Neu03[ IDL09| . 
and not entirely surprising, as it is well known both in models and in field studies of 
semi arid systems that also in uncontrolled systems low rainfall levels lead to patterned 
(or patchy) vegetation, which often can be sustained at lower rainfall levels than a uniform 
state [S/vH Mm 1 rR,DdR,vdKn4l IZKY+1 3) : this is precisely what happens here as well. How¬ 
ever, the quantitative differences between steady states of Q and (19) are quite significant, 
in particular in the sense (c),(d): The controlled system sustains vegetation at much lower R 
values, and at equal R it yields a much higher profit than (19). Moreover, the computation of 
canonical paths from some initial state (often a CSS) to some OSS yields precise information 
how to go to the OSS to maximize the objective function. 


Remark 1.4. In economics, the co-states A and /i are also called shadow prices^ which are 
sometimes difficult to interpret; here, however, the shadow price A has a nice interpretation 
as an optimal tax for private optimization, as follows. Introducing a tax r per unit harvest, 
i.e., setting 


n{v, E) = {p- r)v^E^-^ - cE, 


( 21 ) 


we obtain 


E = 


(p-r)(l -o;)y/" 


( 22 ) 


instead of (20). Thus, after solving the optimization problem via Pontryagin’s Maximum 
Principle, and thus in particular computing the shadow prize A of the vegetation constraint 
along an optimal path, comparing (22) to we see that private optimization maximizes J 
from if the tax r is set in an in general time and space dependent way as r = A. J 


As already said, our method follows |GU15j . where we study optimal controls for a model 
of a shallow lake ecology/economy, given by a scalar parabolic PDE. However, the results are 
rather different. First of all, without control, i.e., for a fixed control equal to some parameter, 
the scalar PDE in |GU15j shows no pattern formation: the patterns in |GU15j are only due to 
the control, which may be called “control induced pattern formation” or, as in |BX08[lBX10j . 
“optimal diffusion instability” (GDI). However, the parameters in |GU15| have to be carefully 
fine-tuned to obtain POSS, which moreover are only locally stable, see also |Gral5| . Here, 
to obtain POSS we need no fine-tuning of parameters. 

From the methodical and algorithmic point of view, our results for ^ illustrate that 
our two-step approach is well suited to deal with non-uniqueness of CSS in nonlinear PDE 
optimal control problems, and the typically associated multiple canonical paths and multiple 
local maxima of the objective function. See also ^for further discussion of efficiency. 
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2 Numerical simulations 


2.1 ID canonical steady states bifurcating from the FCSS branch 

The bifurcation scenario for the stationary problem G{u) = 0 can be studied conveniently 
with pde2path. First we concentrate on the ID case u = u{x)^ x G (—T, L), where the domain 
length must be chosen in such a way to capture pertinent instabilities of the FCSS branch. 
In |BX08|. IBXIO] conditions for pattern forming instabilities in terms of the Hamiltonian T-L 
and its derivatives at a FCSS are given. These are similar to the well known Turing space 
conditions |Mur89| . and moreover allow the calculation of the critical wave-number kc of the 
bifurcation patterns. For instance, at i? = 5 |BX10| find (A:_,A;+) (0.146,1.455) for the 

band of unstable wave numbers at the FCSS. 

If one is interested in accurately capturing the first bifurcation, then one should either fit 
the domain to the (wave number of the) first instability (see, e.g., |UW14j for examples), or 
use a very large domain, which gives a rather dense set of allowed wave numbers. However, 
for simplicity, and with the (expensive) t-dependent canonical paths in mind, here we do 
not want to use a very large domain, and, moreover, rather take the point of view that the 
domain comes with the model. Thus, we do not want to be too precise on fitting the domain 
to the first instability over an infinite domain, and simply choose L = 5. Of course, increasing 
the domain size (certainly in integer multiples of L) will only increase the number of patterns 
and bifurcations, and on the other hand there is a critical minimal domain size below which 
no patterns exist. 

In order to present our results in a domain independent way we give averaged quantities 
such as Jc,a 5 see ( [TT| ), and 

(v) := f v{x) dx, and so on. (23) 

Jn 

Figure [^a) shows a basic bifurcation diagram of ID CSS. We start with a FCSS at i? = 34 
which can be easily found from an initial guess ('u, re. A,/i) = (400,10,0.5,1) followed by a 
Newton loopj^ See Table § for numerical values of some of the FCSS found this way, and of 
other selected points in the bifurcation diagrams. Following the FCSS branch with decreasing 
R we find a number of branch points to PCSS, and near i? = 4 we find a fold for the FCSS 
branch. The lower FCSS branch continues back to large i?, but is not interesting from a 
modeling point of view. The upper FCSS branch has the SPP until the first bifurcation at 
R — Rc 21.5, where a PCSS branch pi with period 20/3 bifurcates subcritically; see the 
example plots of solutions on pi. This is a pitchfork bifurcation, but here and in general we 
only follow the branch in one direction; the other direction is often related to the first by 
symmetry, e.g., spatial translation by half a period. 

The pi branch has a fold at 31, after which it has the SPP down to a secondary 

bifurcation at i ?2 ~ 9.4. Near i? = 3.1 the pi branch also has a second fold, after which 
it continues to larger i? as a branch of small amplitude patterns, (a) also shows the PCSS 

^([^,b) also has the trivial solution branch (u, w) — (0, R^/tw)^ which yields the trivial branch FCSSo with 
(u, w, A, /i) = (0, 0, 0) (and hence E = 0), which however is of little interest here. 
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(a) partial bifurcation diagram of CSS 



(d) example CSS plots, with costates 

p1/pt11,J -22 






(b) Jc,a over R 



R 




(c) example CSS plots 


p2/pt23, Jj, ^-13.61 



- 10w 

- 10E ' 




o'-^- 

-5 0 5 





Figure 1: (a) (partial) bifurcation diagram of CSS in ID, including a (small) selection of bifurcation 
points indicated by o. The labeled points are tabulated in Table and a selection of these solutions 
together with E (and the co-states for the pl-branch) are plotted in (c),(d). (b) shows the two 

branches FCSS and pi with colors as in (a) in a Jc,a over R diagram, which allows to identify which 
CSS maximizes Jc^a amongst the CSS (all other branches from (a) have a significantly lower Jc,a)- 


branches bifurcating from the second and third bifurcation points on the FCSS; interestingly, 
p3 connects back to p2 in a secondary bifurcation. However, except for the FCSS branch for 
R > Rc and the pi branch between the first fold and the first secondary bifurcation, no other 
branch has solutions with the SPP, cf. (16). 

Figure [^b) shows the FCSS and pi branches in a Jc,a over R bifurcation diagram. This 
already shows that at, e.g., R — 20 (in fact, for R smaller than about 24.4), solutions on 
pi yield a larger Jc^a than the FCSS, and thus appear as a candidates for POSS. From the 
applied point of view, probably the most interesting aspect of the solution plots in (c) and 
(d) is that after the first fold on the pi branch the effort E has local minima, not maxima, at 
the maxima of v. Instead, E has maxima on the slopes near the maximum of v. Taking into 
account and the co-state plots in the second row of (d), this can be attributed to the 
distinctive peaks in A at the maxima of v. These peaks in the shadow price of the vegetation 
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point 

FCSS/13 

pl/16 

pl/34 

FCSS/17 pl/11 

pl/38 

FCSS/28 

pl/49 

p2/23 

R 

28 

28 

28 

26 

26 

26 

20 

20 

20 

{v) 

376.32 

337.92 

283.31 

335.36 

311.63 

252.46 

223.59 

175.49 

185.08 

{w) 

9.25 

10.37 

13.17 

9.3 

10.05 

13.34 

9.62 

13.28 

11.64 

(A) 

0.59 

0.53 

0.35 

0.59 

0.56 

0.34 

0.62 

0.34 

0.45 

(M> 

1.09 

1.06 

1.06 

1.04 

1.02 

1.03 

0.87 

0.92 

0.89 

Jc 

25.85 

25.06 

25.02 

22.45 

22 

22.14 

13.66 

14.43 

13.61 

d 

0 

2 

0 

0 

2 

0 

2 

0 

1 

point 

FCSS/44 

FCSS/91 

pl/65 

FSS7121 

PI 797 

p6730 

FCSS/105 

ql/76 

q6/69 

R 

10 

10 

10 

60 

60 

60 

60 

20 

20 

{v) 

75.08 

0.21 

71.33 

79.73 

163.42 

86.5 

1304.5 

183.14 

151.22 

{w) 

11.46 

88.21 

12.71 

65.51 

39.21 

61.4 

10.06 

12.02 

14.54 

(A) 

0.68 

0.9 

0.42 

NA 

NA 

NA 

0.49 

0.35 

0.31 

(m> 

0.5 

0.0006 

0.65 

NA 

NA 

NA 

1.75 

0.91 

0.91 

Jc 

3.51 

0.002 

4.36 

14.3 

29.31 

15.51 

120.8 

13.93 

13.93 

d 

4 

4 

0 

NA(u) 

NA(s) 

NA(u) 

0 

0 

0 


Table 2: Characteristics of selected points marked in Fig. Bid, top half and bottom left in the table), 
Fig. (bottom center in the table, where * denotes values from the case of private optimization) and 
Fig. (2D case, bottom right in the table). NA for the * values means that these values are not 
defined; for the defect the additional u,s are used to indicate unstable vs stable solutions. 


evolution illustrate that it is not optimal to harvest at the peaks in v as this will strongly 
decrease future income. Also note that the (average) shadow prices (A) on the pi branch 
after the fold are lower than on the FCSS branch at the same i?, while at least at low i?, (/i) 
is higher on pi than on FCSS. 

The vegetation patterns (pi branch) survive for lower R (up to i?crit ~ 3.1) than the 
FCSS branch (i?crit ^ 3.7). However, the difference is not large, and this bottom end of R 
will not be our interest here, despite its significance for critical transitions. Instead, we are 
interested in the optimality of CSS for intermediate i?min < i? < i?/, with i?min = 5, say. 


Remark 2.1. Although our picture of CSS obtained above is already somewhat complicated, 
naturally it is far from complete. Firstly, we only followed the first three bifurcations from 
the FCSS branch, and secondly, there are (plenty of) secondary bifurcations on the branches 
pi, p2 and p3, which here we do not follow. In particular, given that the 1.5-modal (in, 
e.g., v) branch pi maximizes J amongst the CSS, a natural question is whether there also 
exist unimodal or 0.5-modal branches, which might given even higher J. The answer is 
(partly) yes: while we could not find a 0.5-modal branch, there is a unimodal branch pO, 
which bifurcates from p2 in a secondary bifurcation, or, more precisely, connects pi and 
p2. See §2.5. 2| for details, where inter alia we study the bifurcation behavior in the price p. 
Moreover, pO then maximizes J amongst the CSS, and, loosely speaking, turns out to be a 
global maximum for Q. 

Nevertheless, until §2.5.2 , for the sake of clarity we restrict to the primary branches which 
bifurcate from the FCSS when varying R. However, one should keep in mind that whatever 
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method one uses to study optimization problems like it is always possible to be stuck 
with some local maxima, and to miss some global maximum. J 


2.2 Comparison to private optimization 

As already explained in the Introduction, private objectives, i.e., individual ranchers maxi¬ 
mizing 7 r('?;, E) = pyOi^i-a _ leads to the system 


dtv = diAv -h {gwv^ — d{l + 6v) — A)v, 
dtw = d2Aw + R{f3 + - {tuV + r^)w, 


(24a) 

(24b) 


with A = 7 ^“^ = 0.543 for the economic parameters (c,p, a) = (1,1.1,0.3) from Table In 
this section we compare the bifurcation diagram in R for steady states of ( [2^ , see Fig. to 
that for Q in Fig|^ 


Roughly speaking, both are similar, but for (24) the bifurcations to patterned steady 


states occur at larger i?, and of course also have to be interpreted differently. First of all, we 
start the bifurcation diagram at i? = 130 with a dynamically stable ffat steady state (FSS) 
of (24), which loses stability at Rc ^ 122 due to a supercritical pitchfork bifurcation of a 


branch pine of patterns with period 5. There are a number of further bifurcations from the 
FSS branch; as an example we give p6nc. The pine solutions lose stability in a secondary 
bifurcation near i? = 61 (not followed here), and eventually all branches undergo a fold 
between i? = 36 and R = 24, and turn into small amplitude branches. 

Similar to the CSS case, here we also have (7r(plnc)) > (7r(FSS)) i.e., the patterned states 
yield a higher (average) profit than the ffat states In (c,d) we compare the FSS branch 
with the FCSS branch. Besides again showing that the fold of the FCSS is at much lower i?, 
and hence the socially controlled system supports a uniform vegetation down to much lower 
i?, this also illustrates that, at given i?, (P) and {v) are significantly higher on the FCSS 


branch. Finally, (24) has the trivial branch (v^w) = (0, P^S/r^^;), which however again is of 
no interest to us. 


2.3 Canonical paths 
2.3.1 Main results 


Having computed CSS branches is only half of the program outlined above; we also need to 
solve the time dependent problem (18) to 

• given a point (vq^wq) and a CSS fi, determine if there exists a canonical path from 
('?;o,ieo) to u; 

• compare canonical paths t u{t) (or {v^w^E){t)) to different i.e.: compute and 


compare their economic values J{u) J{v^^wq^E)^ cf. 10, to find optimal paths. 
Assuming that the spatial mesh consists of n points, we summarize the algorithm for the 
first point as follows, with more details given in |GU15[ IUecl5| . First we compute T G 


^although Jc,a and (tt) have different interpretations, in Tablej^we use Jc,a also for (tt), as both are actually 
defined by the same expression 
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(a) (b) 


\r 

/ 


V 

2w 

V-Z 

2E 

5 0 5 




W 

^ /\\ 


V 

2w 

2E 


(c) (d) 




Figure 2: (a) Bifurcation diagram for the case of private optimization, where (tt) denotes the average 
private profit (given by the same formula as the profit P under control); the blue branch is the primary 
bifurcation pine, the green is p6nc. (b) example solutions from (a), plnc/pt97 (top) and plnc/pt30 
(bottom). (c,d) comparison of the flat states without control (black) with the OC case (blue). 


m2' 


nx4n 


corresponding to the unstable eigenspace of u such that the right BC u{T) G Es{u) 


is equivalent to — u) — 0. Then, to solve (18) we use a modification mtom of the two- 

point BVP solver TOM |MT04j . which allows to handle systems of the form M^u — —G{u)^ 
where M G is the mass matrix arising in the FEM discretization of ([^. A crucial step 

in solving (nonlinear) BVPs is a good initial guess for a solution t i-G u{t)^ and we combine 
mtom with a continuation algorithm in the initial states, again see |GU15[ §2.1] for further 
discussion, and |Uecl5| for implementation details. 

For the second point we note that for a CSS u we simply have J{u) = Jc,a{u)/p. Given a 
canonical path u{t) that converges to a CSS fi, and a final time T, we may then approximate 

. e-p^ 

J{vo,Wo,E)= eP*Jc,a{v{t,-),E{t,-))dt-\ - Jc,a{u). (25) 

Jo P 


In principle, given u — with d{u) = 0 we could choose any ('L’ 0 ,'^ 0 ) in n 

neighborhood of (f), w) (or globally, if is a globally stable OSS) and aim to find a canonical 
paths from ('L’ 0 ,'^ 0 ) lo u. However, the philosophy of our simulations rather is to start at the 
state values of some CSS, and see if we can find canonical paths to some other CSS which 
give a higher J. We discuss such canonical paths in decreasing i?, starting with i? = 26 in 
Fig. and postponing the situation at i? = 28 to §2.3.2[ 


Remark 2.2. a) Note again, that our discussion is based on the primary bifurcations in R 
from the FCSS in Fig. which misses a branch of unimodal CSS, cf. Remark 2T and §2.5.2 


b) Although only the state values (vq^wq) are fixed as initial conditions for canonical paths 
as connecting orbits, in order not to clutter notations and language we write, e.g., i^fcss^ps 
for a connecting orbit starting at the state values of fiFCSS sind connecting to ups- \ 
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R = 26. At i ? = 26 we have two CSS with d{u) = 0, namely ufcss given by FSS/ptl7, 
and tips given by pl/pt38. Figure shows two canonical paths to these CSS. (a) shows 
ll'L’(t) —'L’olloo —'L’llloo for fho path from the “intermediate” patterned CSS tipsi given 

by pl/ptll to tipcss- This indicates that and how fast the canonical path leaves the initial 
{vq^wq) and approaches the goal = (t),re) (the differences in the second component 

\\w{t) — re*||oo are always smaller). Moreover we plot 4Jc,a(^), illustrating that Jc,a{i) does 
not vary much along the canonical path. However, the differences may accumulate. 

(a) Path from upsi to liFCSs: diagnostics, control F, vegetation v and water w. 
p1/pt11 toFSS/pt17 



(b) Path from upsi to hps: diagnostics, control F, vegetation v and water w. 


FSS/pt17 to p1/pt38 



(c) Co-state paths from PSI to FCSS (left), and from FCSS to PS (right) 



Figure 3: R = 26. (a) convergence behavior, the current value profit, and obtained objective values 
for the canonical paths from pl/ptll to the FCSS (left), and F,u,u; on the path; (b) the same for 
the canonical path from the FCSS to pl/pt38. (c) co-state paths from (a),(b). 

The values of the solutions are as follows. We readily have 

J({ipsi) = 733.37 < J{ups) = 738.12 < J(%css) = 748.29, (26) 

and for the paths r^psi^FCSS (^^b), i^fcss^ps (c), and i^psi^ps (i^of shown) we obtain 

T(ttpsi^FCSs) = 735.51, J{upsi^ps) = 744.28, and J(i^fcss^ps) = 749.53. (27) 
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The result for J(i^psi^FCSs) seems natural, as controlling the system to a CSS with a higher 
value should increase J. However, the results for J(i^psi^ps) -^('^fcss^ps) fii'st 

seem counter-intuitive. In i^fcss^ps we go to a CSS with a smaller value, but the transients 
yield a higher profit for the path. In particular, this shows that a CSS which maximizes J 
amongst all CSS is not necessarily optimal Similarly, although ups as a CSS has a lower 
value than upcss^ starting at tipsi it is advantageous to go to ups rather than to liFCSS- 

Due to folds in the continuation in the initial states, again see |GU15| for details, we could 
not compute a path from tips fo tipcss- Thus we conclude that such paths do not exist, and 
(tentatively, see Remark 2 . 2 ) classify tips as an at least locally stable POSS, with tipcss and 
tipsi in its domain of attraction. 

The control to go from tipsi to tipcss (b) is intuitively clear: Increase/decrease E near 
the maxima/minima of vq. Going from tipcss fo tips (c) warrants a bit more discussion: 
For a short transient, E is reduced around the locations X 2 = — 5 and X 4 = 5/3 of the 
maxima of pps- This is enough to give an increase of v around X 2 , 4 . However, under the 
given conditions this does not decrease soil water near ^ 2 , 4 , i.e., the increased infiltration at 
larger v dominates the higher uptake by plants. After this transient, E increases near ^ 2 , 4 , 
thus producing the higher J; see also the discussion of Fig. below. As the behavior of E 


follows from (6 
in (a), (b) in 


, i.e., E 

k He). 


_ AXi a) ^ ^ we also plot A, /i for the paths 


(a) diagnostics 


(b) variables on canonical path 


FCSS/pt44 to p1/pt65 
Jj^=116.9, J=132.49, J^=i45.26 



t 



Figure 4: The canonical path from FCSS to PS at R = 10. 


Smaller R. For i?<i?c~21.5, in Fig|^the only CSS with the SPP is the patterned state 
ups{R)- In Fig. i we focus on the case i?= 10 , and here only remark that the results are 
qualitatively similar for all We show the canonical path from the FCSS to the 

PCSS, where now the strategy to reach a patterned state already indicated in Fi^^c) becomes 
more prominent. Up to t ^ 10, E is reduced around X 2 = —5 and X 4 = 5/3. Conversely, E 
is initially increased near the minima xi = —5/3 and X 3 = 5 of f), leading to a decrease of 

V and an increase of w near ^ 1 ^ 3 . On the other hand, due to diffusion of w the increase of 

V near ^ 2,4 does not lead to a decrease of w compared to the FCSS wq. Instead w increases 
significantly everywhere. After this transient the harvesting effort E is increased near ^ 2 , 4 , 
leading to an overall quick convergence of u{t) to the PCSS u. 
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Thus, the main point for the strategy to go from ufcss fo tips is to initially introduce a 
spatial variation (of the right wavelength) into which yields maxima of v at the minima 
of this initial but then to rather quickly switch to the harvesting on the slopes of the 
generated maxima of v. The canonical path shows precisely how to do this. Also note (blue 
curve in (Fig. [^)) that the initial harvesting briefly yields a higher Jc,a than Jc,a(rips) but 
in the transition 4 < t < 25, say, Jc,a(t) is signiflcantly below Jc,a(fips)- 
For the values we have 


J{ufcss) = 116.9 < J(rtFcss^Ps) = 132.49 < J(rips) = 145.26. 


(28) 


Thus, again tentatively, see Remark |2.2| and in particular (32) below, we classify rips sit 
i? = 10 as a POSS, with ripcss iri its domain of attraction. For completeness we remark that 
at i? = 20 we have 


-^(ripcss) = 455.31, J(rips) = 480.88, and J{ufcss^ps) = 474.57. (29) 


2.3.2 A patterned Skiba point 

In ODE OC applications, if there are several locally stable OSS, then often an important 
issue is to identify their domains of attractions. These are separated by so called threshold 
or Skiba-points (if A = 1) or Skiba-manifolds (if N > 1) |Ski78[ lOCF^OSi [KW10| . Roughly 


(a) A Skiba at a 


0.9 (b) Paths of (almost) equal values to the FCSS and the upper PCSS. 



Figure 5: In (a), the blue and red lines gives J for the canonical paths u^fcss sind u^ps from 
{v^w)a{0) := a{v^w)ps + (1 — (^){v^w)pcss to FCSS/ptl3 and to pl/pt34, respectively. In (b), the 
white surfaces are for u^pcss rind the colored ones for rt^ps- 


speaking, these are (consist of) initial states from which there are more than one optimal 
paths with the same value, but eonneeting to different CSS. In PDE applications, even under 
spatial discretization with moderate nN, Skiba manifolds should be expected to become very 
complicated objects. 

In Fig. we just given one example of a patterned Skiba point “between” ups given by 
pl/pt34 and upss given by FCSS/ptl3, at i? = 28, where tips ^^nd tipcss ^^re the two possible 
targets for canonical paths. The blue and red lines in (a) gives J for the canonical paths as 
indicated. The lines intersect near a ^ 0.9, giving the same value J. Hence, while the two 
paths are completely different, they both are equally optimal, and for illustration (b) shows 
the two paths for a — 0.9, where | J-^fcss ~ -^-^ps| < 0.08. 
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2.4 2D results 


The basic mechanisms of pattern formation in reaction-diffusion models can usually be stud¬ 
ied in ID, but for quantitative results for vegetation-water ecosystem models one should also 
consider the more pertinent 2D situation, and clearly the same applies to the OC system. 
Even though we do not expect qualitatively different results from the ID case, here we give 
a short overview over 2D PCSS and the associated canonical paths. 


(a) bifurcation diagram (b) Jc,a over R 




(c) example CSS 

V at q1/pt76 



E at q1/pt76 

-5 0 


V at q6/pt69 



-5 0 5 



Figure 6: Partial bifurcation diagram and example plots of CSS in 2D. 

The first question concerns the second spatial length, now called ^-direction. For classical 
Turing bifurcations, typically stripes and/or hexagonal spots are the most stable bifurcating 
2D patterns, see |UW14| and the references therein for discussion, and thus by analogy here 
we choose the domain Q = (—L,L) x ( —v^L/2, v^L/2) with L = 5 as in ID. Figure]^ a) 
shows five branches bifurcating from the FCSS. It turns out, that the ID branch pi actually 
comes out of the 5th bifurcation point in 2D, and is therefore called q5 now, while the first 
branch ql corresponds to horizontal stripes, see (c). Thus, L = 5, chosen for simplicity, does 
not capture the first instability in ID. However, the first bifurcation points are very close 
together; moreover, upon continuation to low R the q5 branch still yields the highest 
see (b). 

The sixth branch q6 is a hexagon branch. Similar to q5, both ql and q6 have the SPP after 
their first folds. The other branches are other types of stripes or spots, for instance “squares”, 
but none of these have the SPP. All branches exhibit some secondary bifurcations, not shown 
here, and to not overload the bifurcation diagram we only plot the starting segments of q2 
(green) and q3 (magenta). 

At, for instance, i? = 20 we now have 3 possible targets for canonical paths: flhs (horizontal 
stripes) from ql, lihex (hexagons) from q6, and (vertical stripes) from q5, already discussed 
in ID as pi. It turns out that we can reach each of these from the FCSS, with J('Ufcss) = 
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t=3.1 



Figure 7: Snapshots of E and (v^w) on the canonical path from the FCSS to the hexagons. The 
states v^w both go directly into their ultimate hex pattern, which then only grows, while E shows a 
switching analogously to the ID case. 


455.31, cf. (29), and similarly J{uvs) — 480.88 and liFCSS^vs = 474.56 as in ID (since these 
values are normalized by |D|). For the path to lihs with J{uhs) — 464.33 we obtain J — 464.91, 
and the path to fihex with J(iihex) = 464.33 = J(iihs) (up to 2 digits) yields J — 467.38. Thus 
we again have V{ufcss) = 474.6. The strategies for these paths are the natural extensions 
of the ID case: given a target PCSS initially E has minima at the maxima of i), but after 
a rather short transient during which •) develops maxima at the right places, E changes 
to harvesting in the neighborhood of these maxima. Movies of these paths can be found at 
and in Fig. we present some snapshots. We could not compute canonical paths 
from any of the PCSS to any other PCSS, with the continuation typically failing due to a 
sequence of folds. Thus we strongly expect all three PCSS to be locally stable POSS. 


2.5 Remarks on further parameter dependence 

So far we varied the rainfall R as our external bifurcation parameter. Similarly, we could 
vary some other of the physical parameters 77 ,..., di^2 (first six rows of Table [^, and in 
most cases may expect bifurcations to patterned CSS. 

Maybe even more interesting from an application point of view is the dependence on the 
economic parameters p, c,p and a (discount rate, cost for harvesting, price of harvest, and 
elasticity), as these may vary strongly with economic circumstances. Moreover, varying a 
second parameter often also gives bifurcations to branches which were missed upon contin- 


20 







nation of just the primary parameter, and these may play an important role in the overall 
structure of the solution set; this does happen here, see §2.5.2 below. 


2.5.1 Experiments with the discount rate p 

In Figj^we illustrate the dependence of the PCSS on the pi branch from Fig. [T]on p, at fixed 
R = 10, cf. also Fig. Panel (a) shows the bifurcation diagram; to obtain the blue and black 
branches we continued the points pl/pt65 and FCSS/pt44 from p = 0.03 down to p = 0.005, 
reset counters, renamed pi to rl, and continued to larger p again. Both, the FCSS and the rl 
branches then show some folds at p 0.185 and p 0.325, respectively. More importantly, 
the rl branch has the SPP for small p, but loses it at p 0.032 to another PCSS branch si. 
Solutions on si have maxima of different heights in see (b), and have the SPP only up to 
the fold at p = pj ~ 0.046. Moreover, there are further bifurcations from the FCSS branch 
to PCSS branches, but none of these has the SPP. 


(a) Bifurcation diagram in p (b) example plots on rl and si branches 



Figure 8: Bifurcations when varying the discount rate p at = 10. The blue branch (which at 
p = 0.03 corresponds to Fig. loses the SPP at p ~ 0.032, where the red sl-branch bifurcates, which 
has the SPP up to its fold. 

With ups given by sl/pt9 (p = 0.04), the control to go from the FCSS to u shows a 
similar switching strategy as, e.g., in Fig. The values are 

-^('^Fcss) = 82.5, J{upcss^ps) = 92.8, J{ups) = 106.8, (30) 

which are more or less comparable to Fig. (with p = 0.03). On the other hand, for the 
canonical path from the FCSS to rl/ptO at p = 0.005 we obtain 

J{upcss) = 774.9, J{upcss^ps) = 892.2, J{ups) = 897.7. (31) 

Additional to the larger total values due to the smaller discount rate, compared to Fig. the 
canonical path to ups now has almost the same value as tips itself. This illustrates that at 
low p the transients have less influence, as expected. 
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For p > pf none of the CSS on the branches that are shown in Fig. or that can be 
obtained from the shown bifurcation points, have the SPP. This does not mean that PCSS 
with the SPP do not exist in this parameter regime, but rather that they must be obtained 
by continuation and bifurcation in some other way, cf. Remark |2.1| and §2.5.2 

2.5.2 Dependence on the price p, and the unimodal branch 

In Fig. a) we illustrate the dependence of the FCSS and pi branches on the price p, with 
fixed R = 10, starting at p = 1.1 with pl/pt65 and FCSS/pt44 from Fig. respectively. 
Naturally, the values decrease as p decreases, and not surprisingly pi bifurcates from the FSS 
branch at some Pc~0.55. Next, as an additional benefit we find the “unimodal” branch pO, 
which bifurcates from the FCSS branch near p=0.5, and which yields a higher Jc^a than pi. 


(a) Jc^a over p (b) Continuation of pO in R 




(c) Example plots on pO and pOl 




p01/pt38,J^ ^=14.05 

6001- /-v 
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Figure 9: (a) branches FCSS (black), pi (blue), and pO (red) over p, R = 10 fixed, (b) Continuation 
in R of new branch pO found in (a), and bifurcating branch pOl (green), together with known branches 
FSS (black) and p2 from Fig. Example plots in (c). 


Thus, in a next step we continue pO to p = 1.1 and then switch back to continuation in i?, 
with p = 1.1 fixed, i.e., pO/ptO in (b)-(c) is ptll2 from (a). It turns out that the pO branch 
has the SPP up the fold at i?f 30.15, and slightly below the fold there is a bifurcation 
point to the green branch. This contains some “skewed” solutions, and connects pO and pi. 
Ultimately, pO connects back to p2 from Fig. [T] at low amplitude near R = 21.1. Thus, we 
could also have found pO by following secondary bifurcations in Fig. The values pertinent 
for the canonical path from FCSS/pt44 to pO/ptO are 

-^(f^FCSs) = 116.9 (as in ([^), T(ttFCSS^Po) = 145.11, J(fipo) = 165.42, (32) 

which shows that the path to pO dominates the path to pi from Fig. Thus, the point 
FCSS/pt44 is in the domain of attraction of pO, and not of pi. 

On the other hand, we could not find canonical paths from pl/pt65 to pO/ptO (or vice 
versa). Therefore pl/pt65 can still be classified as an at least locally stable POSS, and 
similarly, pO/ptO is only locally stable^ since it does not attract pl/pt65. Next one could 
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compute a number of Skiba-points (cf. p.3.2 ) “between” pO and pi to roughly characterize 
the respective domains of attraction, but here we refrain from this. 

Finally, despite trying some further combinations of continuations/bifurcations and also 
suitable direct initial guesses followed by a Newton loop, we could not find a “half”-modal 
branch in the parameter regions studied so far, i.e., a branch on which solutions are monotonous 
in V. Thus, it appears that such “long wave” PCSS, i.e., with period 20, do not exist in these 
parameter regions. 


3 Summary and discussion 

Our numerical approach for spatially distributed optimal control problems may yield rich 
results, if applied carefully, in the following sense. First, the canonical system may have 
many steady states, and it is in general not clear how to find all relevant CSS, and which 
of the CSS have the SPP and hence are suitable targets for canonical paths, and second it 
needs to be checked which CSS ultimately belong to optimal paths. 

On the other hand, the value J of the CSS can easily be calculated in parallel with the 
bifurcation diagram of the CSS, allowing to identify which CSS maximize J amongst all CSS. 
Compared to the computation of canonical paths (or direct methods for the optimization 
problem Q), this first step is relatively cheap numerically, but (together with the SPP) 
typically still gives a strong indication for optimal CSS. 

The computation of canonical paths is a connecting orbits problem, and in particular in 
two spatial dimensions this may become numerically expensive. In practice we found our 
two-step approach to be reasonably fast for up to 5000 degrees of freedom of u at fixed time, 
e.g., for our vegetation model 1250 spatial discretization points and 4 components, and up to 
100 temporal discretization points; up to such values a continuation step in the calculation of 
a canonical path takes up to a minute on a desktop computer (Intel i7, 2.3 GHz), such that a 
typical canonical path is computed in about 5 continuation steps in at most 5 minutes, and 
often much more quickly. In ID, with n — 50, say, typical canonical paths are computed in 
a few seconds. 

Here we applied our method exemplarily to the optimal control model from |BX10j . and 
for this we again summarize our main results as follows, cf. also (a)—(d) in the Introduction: 

(i) Compared to the case of private optimization, we have CSS (both flat and patterned) for 
significantly lower rainfall values i?, and the whole Turing (like) bifurcation scenario 
is shifted to lower R. This is important for welfare as it means a much increased 
robustness of vegetation (and hence harvest) with respect to low rainfall. Moreover, at 
a given R the social control gives a significantly (almost ten times) higher J for steady 
states than private optimization; see, e.g.. Table bottom center, and Fig. He), for 
numerical values. 

(ii) At low i?, some PCSS yield a higher J than the FCSS, and some of these PCSS are 
locally stable POSS. 

(hi) The optimal controls to reach such a POSS u from a FCSS follow some general rules: 
first decrease the harvesting effort E at the location of the maxima of the desired POSS, 
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then start harvesting near (but not at) the maxima of the POSS, as determined by E 
from the POSS. The increase in welfare by controlling the system from a FCSS to a 
POSS can be up to 40%, see (32). See also, e.g., (27), (28), (29), and the values given 
in {2^ for further numerical values. 

(iv) The co-state (or shadow price) A computed for an optimal path can be interpreted as 
an optimal tax for private optimization, see Remark |1.4[ 

The points (ii) and (iii) emphasize that resource management rules in systems with multiple 
CSS may need to take several of these into account, as also illustrated by the Skiba point 
computations in §2.3.2 We strongly expect similar results about POSS in other spatially 
distributed optimal control problems for (Turing like) systems of PDE. Thus we hope that 
our numerical approach is a valuable tool to study the basic behavior of spatially distributed 
optimal control models with the states fulfilling a reaction diffusion system. 
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